Locating periodic orbits by Topological Degree theory 
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Q ■ Abstract 

We consider methods based on the topological degree theory to compute periodic 
orbits of area preserving maps. Numerical approximations of the Kronecker integral 
and the application of Stenger's method allows us to compute the value of the 
topological degree in a bounded region of the phase space. If the topological degree 
of an appropriate set of equations has a non-zero value, we know that there exists 
at least one periodic orbit of a given period in the given region. We discuss in 
detail the problems that these methods face, due to the existence of periodic orbits 
near the domain's boundary and due to the discontinuity curves that appear in 
maps defined on the torus. We use the characteristic bisection method for actually 
locating periodic orbits. We apply this method successfully, both to the standard 
map, which is a map defined on the torus, and to the beam-beam map which is a 
continuous map on the plane. Specifically we find a large number of periodic orbits 
of periods up to 40, which give us a clear picture of the dynamics of both maps. 
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> ! 1 The topological degree (TD) and its computation 
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We consider the problem of finding the solutions of a system of nonlinear equations 
of the form 

F n (x) = e n , (i) 

where F n = (fi, f 2 , . . . , f n ) '■ D n C R n K n is a function from a domain D n into 
R n , 6 n = (0, 0, . . . , 0) and x = (xi, x 2 , . . . , x n ). The above system can be written 
as: 

f 1 (x ll x 2 , . . . ,x n ) = 0, 
f 2 (x u x 2 ,...,x n ) = 0, 

f 3 (x 1 ,x 2 ,...,x n ) = 0. 

The topological degree (TD) theory gives us information on the existence of solu- 
tions of the above system, their number and their nature || |3], |(|. Kronecker 
introduced the concept of the TD in 1869 ||, while Picard in 1892 || succeeded 
in providing a theorem for computing the exact number of solutions of system 
(H). Numerical methods based on the TD theory have been applied successfully to 
numerous dynamical systems (e.g. |16|, [17|, [18], [19], §). 
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In order to define the concept of the topological degree we consider the function 
F n of system (U) to be continuous on the closure D n of D n , satisfying also F n (x) ^ 
Q n for x on the boundary b(D n ) of D n . We also consider the solutions of (jl[) to 
be simple i.e. the determinant of the corresponding Jacobian matrix (Jf„) at the 
solution, to be different from zero. Then the topological degree (TD) of F n at Q n 
relative to D n is defined as: 



deg[F n , D n , n ] 



Yl sgn(detJ Fn (x)) 



(3) 



where detJ^ n is the determinant of the Jacobian matrix of F n , sgn is the well- 
known sign function, iV + the number of roots with det Jp n > and AL the number 
of roots with det Jf„ < 0. It is evident that if a nonzero value of deg[F n , D n , Q n ] 
is obtained then there exist at least one solution of system F n (x) = n within D n 

A practical way to find the TD is the computation of Kronecker integral 
In particular, under the assumptions of the above-mentioned definition of the 
the deg[F n , D n , n ] can be computed as: 



deg[F n , D n , 9 r 
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and T(x) is the gamma function. 

In order to find the number of solutions of system (|l|) we consider the function 



where 
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7 n+l — (fl, f2, ■ ■ ■ , fn, fn+l) '■ Ai+1 C 

f n+1 = y detJ Fn , 
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: xi,X2, ■ ■ ■ ,x n ,y and D n+ i is the product of D n with a real interval on the 
l/-axis containing y = 0. Then the exact number A^ of the solutions of equation 
F n (x) = 9 n is proven to be [H: 



N = deg[F n+1 ,D n+1 ,e 



n+l 



By applying this result and using the computation of TD by Kronecker integral 
) in the case of a set of two equations: 



fi(x 1 ,x 2 ) 
f2(x 1 ,x 2 ) 



0. 
0. 



(9) 



we find that the number A^ of roots in the domain D 2 = [a, b] x [c, d] is given by: 

Qdx\dx 2 
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where e is an arbitrary positive value, and 
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Figure 1: (a) The polyhedron ABDC is noncharacteristic while the polyhedron AEDC is charac- 
teristic, (b) Application of the characteristic bisection method to the characteristic polyhedron 
AEDC, giving rise to the polyhedra GEDC and HEDC, which are also characteristic. 



with J denoting the determinant of the Jacobian matrix of F 2 = (fx, / 2 ). 

Another method for computing the TD of F n at a domain D n is the application 
of Stenger's theorem |12, |7j. Following this approach for finding the TD, we only 
need to know the signs offunctions fx, / 2 , ...,/„ in a 'sufficient' set of points on 
the boundary b(D n ) of D n . 



2 The characteristic bisection method 

The characteristic bisection method is based on the characteristic polyhedron con- 
cept for the computation of roots of the equation ([!]). The construction of a suitable 
n-polyhedron, called the characteristic polyhedron, can be done as follows. Let M n 
be the 2 n xn matrix whose rows are formed by all possible combinations of —1 and 
1. Consider now an oriented n-polyhedron IP, with vertices Vj., k = 1, . . . , 2 n . If 
the 2 n x n matrix of signs associated with F and IP, whose entries are the vectors 

sgnF n (V k ) = ( Sg nfx(V k ),sgnf 2 (V k ),...,sgnf n (V k )), (12) 

is identical to M n , possibly after some permutations of these rows, then IP is called 
the characteristic polyhedron relative to F n . If F n is continuous, then, after some 
suitable assumptions on the boundary of IP we have: 

deg[F n ,IP,e n ] = ±1^0. (13) 

This means that there is at least one solution of system F n (x) = n within IP. 

To clarify the characteristic polyhedron concept we consider a function F 2 = 
(/ij/2)- Each function % = 1,2, separates the space into a number of different 
regions, according to its sign, for some regions fa < and for the rest fi > 0, 
% = 1,2. Thus, in figure |l](a) we distinguish between the regions where f\ < 
and / 2 < 0, fx < and / 2 > 0, fx > and / 2 > 0, fx > and / 2 < 0. Clearly, 
the following combinations of signs are possible: (— , — ), (— , +), (+, +) and (+, — ). 
Picking a point, close to the solution, from each region we construct a characteristic 
polyhedron. In this figure we can perceive a characteristic and a noncharacteristic 
polyhedron IT 2 . For a polyhedron II 2 to be characteristic all the above combinations 
of signs must appear at its vertices. Based on this criterion, polyhedron ABDC 
is not a characteristic polyhedron, whereas AEDC is. A characteristic polyhedron 
can be considered as a translation of Poincare-Miranda hypercube ]T3[ . 



Next, we describe the characteristic bisection method. This method simply 
amounts to constructing another refined characteristic polyhedron, by bisecting a 
known one, say IT n , in order to determine the solution with the desired accuracy. 
We compute the midpoint M of an one- dimensional edge of IP, e.g. (Vi, Vj). The 
endpoints of this one-dimensional line segment are vertices of IP, for which the 
corresponding coordinates of the vectors, sgnF n (Vi) and sgnF n (Vj) differ from each 
other only in one entry. To obtain another characteristic polyhedron we compare 
the sign of F n (M) with that of F n (Vi) and F n (Vj) and substitute M for that vertex 
for which the signs are identical. Subsequently, we reapply the aforementioned 
technique to a different edge (for details we refer the reader to [O, [H]. ||). 

To fully comprehend the characteristic bisection method we illustrate in figure 
|l](b) its repetitive operation on a characteristic polyhedron II 2 . Starting from the 
edge AE we find its midpoint G and then calculate its vector of signs, which is 
(—1,-1). Thus, vertex G replaces A and the new refined polyhedron GEDC, is 
also characteristic. Applying the same procedure, we further refine the polyhedron 
by considering the midpoint H of GC and checking the vector of signs at this point. 
In this case, its vector of signs is (—1, —1), so that vertex G can be replaced by 
vertex H. Consequently, the new refined polyhedron HEDC is also characteristic. 
This procedure continues up to the point that the midpoint of the longest diagonal 
of the refined polyhedron approximates the root within a predetermined accuracy. 

3 Applications 

We apply methods based on the topological degree theory to compute periodic 
orbits of two area preserving maps, the Standard map (SM) M, which is a map 
defined on the torus: 

y- : y-l~s^ X) ' n » d < 1 >. W5), (14) 



and the beam-beam map (BM) JT|, [TlJ| , which is a map defined on 

x' = x cos(27ro;) + {y + 1 — e~ x2 ) sin(27ra;), 
y' = — x sin(27ru;) + (y + 1 — e~ x ) cos(27ru;). 



P 2 . 



(15) 



Given a dynamical map M : {x' = gi(x,y),y' = g2(x,y)}, the periodic points 
of period p are fixed points of the p-iteration M p of the map and the zeroes of the 
function: 

F = M p — I = ( = 9i( x ^y) ~ x Qg\ 

I h = g&fay) -y ' 

where I is the identity. 

One can use a color map to inspect the geometry of function F ([16]) and to locate 
its zeroes. The color map is created by choosing a lattice ofmxm points and 
by associating to each point a color chosen according to the signs of the functions 
fii f'2 as shown in figure ^| A simple algorithm allows to detect the cells, formed 
by the lattice of m x m points, whose vertices have different colors. A cell is a 
candidate to have a zero in its interior if the corresponding topological degree is 
found to be different from zero. In figures |3] and ^ we construct the color map and 
apply the above mentioned algorithm for locating periodic orbits of period 3 for 
the SM and of period 5 for the BM, respectively. In both figures the gray circles 
at the right panels denote the positions of the found periodic orbits. We see that 
for both maps some periodic orbits were not found because some of the four color 
domains close to the fixed point were very thin. On the other hand, due to the 
discontinuity of function F ([TBI), some zeros that do not correspond to real periodic 
orbits were found for the SM~(right panel of figure |3|). 




Figure 2: Sketch of the domains where functions /i and of system (16) have a definite sign 
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Figure 3: Standard map ([lj]) for k — 0.9: color map for p — 3 iterations of the map computed 
on a square of m x m points for m = 512 (left panel); phase plot of the map (right panel). The 
gray circles denote the positions of the zeros of the corresponding function (fiff) . 
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Figure 4: Beam-beam map ( jig ) for ui = 0.21: color map for p = 5 iterations of the map computed 
on a square of m x m pointsfor rn — 512 (left panel); phase plot of the map (right panel). The 
gray circles denote the positions of the zeros of the corresponding function fll6|) . 
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Figure 5: The discontinuity curves of the standard map M ( |l4| ) divide the phase space in five 
continuous regions (I, II, III, IV, V). In each region the computation of the TD can be performed 
accurately. 



For maps defined on the torus like the SM flllD , the computation of the TD 
using Stenger's method or the Kronecker integral fllPP faces problems due to the 
presence of discontinuity curves. Indeed Kronecker integral is defined on a domain 
where the function F fll6|) is continuous. 

For the SM the discontinuity curves are the lines x = —0.5 and y = —0.5. By 
applying the SM map M once these lines are mapped on the curves seen in the right 
panel of figure 0. On the initial phase space there exist also the discontinuity curves 
that will be mapped after one iteration to the lines x = —0.5 and y = —0.5. These 
curves are also plotted in the left panel of figure 0. These curves can be produced 
by applying the inverse SM to the discontinuity lines x = —0.5 and y = —0.5. So 
the discontinuity curves divide the initial phase space in five continuous regions 
marked as I, II, III, IV and V in figure 0. In each region the computation of the 
TD can be performed accurately by Stenger's method or by evaluating Kronecker 
integral. If, however, the boundary of the domain where these procedures are 
applied, cross a discontinuity curve the results we get are not correct (figure 0). 

In order to study the dependence of the procedure for finding the TD in a region 
D, with respect to the distance of a root from the boundary of D, we consider the 
simple map 

f fi(x,y) = y-f + x 
1 f2(x,y) = y 



(A,/ 2 



(17) 



The lines f\ = 0, f 2 = are plotted in figure 0(a). The system F* = (0,0) 
has three roots. The determinant of the corresponding Jacobian matrix (detJj?*) 
is positive for root (0,0) and negative for roots (— \/3, 0) and (y/3,0). We also 

consider a rectangular of the form [—a, 2] x [—2, 2] with a > v^3, shown in figure 
0(a). Since this domain contains the three roots of the system the value of the TD 

is —1. We set a = v3 + e with e > so that the boundary approaches the root 
as e — > 0, as shown by the arrow in figure 0(a). We compute the TD for different 
values of e applying Stenger's method, by using the same number of points m on 
every side of the rectangle. We denote by n gp = 4m the smallest number of grid 
points needed to compute the TD with certainty. In figure 0(b) we plot in log-log 
scale, n gp with respect to e (dashed line). The slope of the curve is almost —1 so 
that m oc e _1 . The same result holds for any map when the boundary approaches 
a root (the solid line in figure 0(b) is obtained for a similar example for the SM 

(0))- 

Despite the problems caused by the discontinuity curves or by roots located 
very close to the domain's boundary, the use of the characteristic bisection method 
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Figure 6: (a) Number of period 1 fixed points Af\ evaluated for the SM (Jl4|) with k = 0.9 using 
the Kronecker integral (|lO|), in a rectangular domain whose upper-side moves, as a function of 
the y coordinate of this side. The rectangle and the discontinuity lines are shown in (b). For 
the various rectangles, M\ should be equal to 1 since they contain only 1 fixed point of period 1, 
namely point (0,0). The two points marked by arrows in (a) where Afi deviates from the correct 
value N\ — 1, correspond to y ~ 0.358 and y « 0.466 respectively, where the upper-side of the 
rectangular crosses the two discontinuity curves in (b). 




Figure 7: (a) Plot of the curves /i = y — ^- +x = 0, fa = y = (b) Dependence of the number of 
grid points n gp , needed for computing the correct value of the TD in a domain, on the distance e 
of a root from the boundary of the domain, for the set of equations of (a) (dashed line) and the 
SM (continuous line). 



Figure 8: Periodic orbits up to period p — 40 for the SM ( |l4| ) for k = 0.9 (left panel) and for 
the BM (£51) for w = 0.14 (right panel). Different gray-scales correspond to periodic orbits with 
different land of stability. 

can locate a big fraction of the real periodic orbits. Actually by applying the 
characteristic bisection method on the cells of a lattice formed by 2000 x 2000 
grid points we were able to compute a sufficient number of the periodic orbits 
with period up to 40, for the SM (figure |S], left panel) and the BM (figure ^ right 
panel). The computed periodic orbits give us a clear picture of the dynamics of 
these maps. 

4 Synopsis 

We have studied the applicability of various numerical methods, based on the topo- 
logical degree theory, for locating high period periodic orbits of 2D area preserving 
maps. 

In particular we have used the Kronecker integral and applied the Stenger's 
method for finding the TD in a bounded region of the phase space. If the TD has 
a non-zero value we know that there exist at least one periodic orbit in the corre- 
sponding region. The computation of the TD for an appropriate set of equations 
allows us to find the exact number of periodic orbits. We also applied the char- 
acteristic bisection method on a mesh in the phase space for locating the various 
fixed points. 

The main advantage of all these methods is that they are not affected by accu- 
racy problems in computing the exact values of the various functions used, as, the 
only computable information needed is the algebraic signs of these values. 

We have applied the above-mentioned methods to 2D symplectic maps defined 
on R 2 and on the torus. The methods for computing the TD are applied to con- 
tinuous regions of the phase space, so their use for maps on the torus is limited to 
regions where no discontinuity curves exist. On the other hand the characteristic 
bisection method proved to be very efficient for all different types of maps, as, it 
allowed us to compute a big amount of the real fixed points of period up to 40 in 
reasonable computational times. 
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